Genetic diversity and population structure of critically endangered Dactylorhiza hatagirea (D. Don) Soo from North-Western Himalayas and implications for conservation

Dactylorhiza hatagirea (D. Don) Soo is medicinally important herb, which is widely used in ayurveda, unani, and folk/traditional medicine system to cure diseases. Due to its immense ethno-botanical properties, the trade of D. hatagirea is estimated to be USD 1 billion/year in India. Unfortunately, due to overexploitation of the herb from the wild, has resulted in dwindling of its populations in their natural habitats, which has led to its critically endangered status. Molecular genetic studies are still scarce in D. hatagirea, therefore, in current study, genetic diversity and population structure analysis was carried out of 10 populations (48 individuals) collected from three cold desert regions (2527 m–3533 m amsl) of Himachal Pradesh. Mean observed heterozygosity (Ho) and expected heterozygosity (He) was recorded 0.185 and 0.158. The maximum values for Fst (fixation index) and Nm (gene flow) were recorded 0.945 at locus KSSR14 and 1.547 at locus KSSR 4 respectively. Mean genetic differentiation (Fst) coefficient was estimated to 0.542. Overall, low levels of genetic diversity was recorded in the populations of D. hatagirea, might be due to habitat specificity (alpine meadows ecosystem; humid laden undulating habitat), restricted distribution and high anthropogenic activities. However, two populations viz., Bathad and Rangrik were recorded with high diversity and largest number of private alleles, stipulates that these populations might have high evolutionary significance and response to selection. Dendrogram analysis revealed that the populations of D. hatagirea were clustered into four major clusters, which was supported by Bayesian based STRUCTURE predictions. Clustering pattern of majority individuals of different populations revealed consistency with their geographic origin. Outcomes of current study reveals the status of genetic diversity and population structure of endangered D. hatagirea, which can be futuristically utilised for appropriate planning of conservation strategies.

the Himalayas is considered one of the critical factor responsible for decline in the availability of medicinally important plants and their populations.
Dactylorhiza hatagirea (D. Don) Soo belongs to Orchidaceae family, is native to India, Pakistan, Bhutan, Afghanistan, Nepal, and Tibet, is of high value due to its medicinal and ornamental properties. The plant is found at an altitudinal range of 2500-5000 masl in alpine meadows and grassy slopes 5 . It is a monocotyledonous, terrestrial, perennial herb with erect hollow stem and bears purple flowers with green bracts. It is a slow growing, habitat specific species and regenerate poorly due to its requirement of association with mycorrhizae fungi for germination 6 . Bioactive compounds of D. hatagirea, which are of medicinal importance are dactyloses A, B, and dactylorhins A to E. Bulbous roots of the plant is used to prepare "Salep" which is utilized to cure wounds, diarrhoea, dysentery, cough, fever, burns and fractures 7 . The plant has high annual demand (5000 tonnes) and price of dried rhizome is ~ 72 US $ per kg, which has resulted in its overexploitation from its natural/wild habitat 4 . Moreover, it has also been observed that D. hatagirea in wild shares habitat with highly traded species like Ophiocordyceps sinensis and with expanding species like Polygonum polystachyum, which is also affecting negatively the number of individuals in a population of D. hatagirea 8 . Biodiversity loss is among the most disastrous environmental problem and world is now facing the challenges required to reduce the rate of loss of biodiversity 9 . Moreover, due to the increase in the demand of natural compound in the healthcare system, natural populations of medicinally important plants like D. hatagirea are under tremendous pressure due to overexploitation, indiscriminate extraction and unsustainable management.
Understanding of population structure and genetic diversity pattern is important for sustainable utilisation and planning of conservation strategies for critically endangered plant species. Moreover, maintaining genetic variation in a population is one of the most important objective for conserving critically threatened and endangered species. Whole genome sequencing of critically endangered species to create genomic resources are paramount to understand the genome-wide diversity between and within species, which is beneficial for conservation practitioners. Unfortunately, planning of conservation strategies in critically endangered D. hatagirea is very challenging due to the non-availability of genetic resources. Till date, only 17 nucleotide sequences 10 and one transcriptome study has been reported for D. hatagirea 11 . Few genetic diversity studies using RAPD and ISSR markers from Ladakh region recorded moderate level of diversity among natural populations 12,13 . Variability among natural populations of D. hatagirea from Uttarakhand was also analysed using morphological, biochemical and isoenzymes analysis 14 . However, no population structure and genetic diversity study has been reported from the cold desert regions of Kullu, Kinnaur, Lahaul and Spiti valley of Himachal Pradesh, where D. hatagirea is among the most exploited medicinally important plant for domestic usage and is under increasing demand for export 15,16 .
In current study, attempts were made to perform extensive geographic sampling of natural populations of D. hatagirea, and samples were collected from 10 locations representing three regions (Kullu, Kinnaur and Lahaul & Spiti) of Himachal Pradesh, India. Fifteen stable and highly polymorphic SSR (Simple Sequence Repeats) markers developed in our previous study were utilised for genetic diversity characterization and population structure analysis of collected natural populations of D. hatagirea 10 . SSRs were utilised as they are the most ubiquitous and hypervariable type of repeat elements, which are abundant in the plant genomes. Moreover, SSR markers have added advantages for precise estimate of genetic diversity over other markers because of their codominant nature, high reproducibility and transferability 17 . Accordingly, realising the medicinal, ornamental and industrial applications of D. hatagirea along with its critically endangered status in the Himalayas, current study was performed with an aim to understand the diversity pattern and population structure of natural populations of D. hatagirea for their sustainable management, effective conservation strategy planning and genetic resource management for future breeding programmes.
Population-wise genetic diversity. Population level genetic diversity parameters were found considerably different across 10 populations of D. hatagirea ( Table 2). The average of percentage of polymorphic loci (PPL) across 10 populations was found moderate (42.67%) which ranged from 20.00% (Chhitkul) to 66.67% (Mane) with only three population revealing PPL more than 50% (Mane, Rangrik and Sichling).  www.nature.com/scientificreports/ between Chhitkul and Bathad to maximum of 3.071 between Shakoli and Chhitkul (Table 3). Furthermore, pair-wise comparative analysis of Nei's genetic diversity values revealed variation of 0.019 between Sichling and Mane and 0.891 between Chhitkul and Bathad with majority of the values between populations lies in the range of 0.1-0.7 on the basis of SSR markers (Table 4). Further, Mantel test analysis stipulated weak positive correlation between genetic and geographic distance (Rxy = 0.347, P = 0.04) indicating isolation-by-distance phenomena in D. hatagirea populations (Fig. 1a). Further, AMOVA analysis revealed that 73% of molecular variance was recorded among populations and 27% individual differentiation within populations (Table 5).
Cluster and population genetic structure analysis. Dendrogram analysis using Neighbor-Joining based hierarchical clustering revealed that the 10 populations of D. hatagirea were clustered into four major clusters i.e. cluster I, II, III and IV. Cluster I and IV represented populations from Spiti valley viz., Mane, Sichling, Giyu in cluster I and Rangrik and Giyu in cluster IV. Interestingly, one individual of Giyu (Spiti valley) population was found intermixed with these two populations. Cluster II represented individuals from Bathad population (Kullu) and cluster III represented majority of individuals from Shakoli (Lahaul valley), Tosh (Kullu), Sangla and Chittkul (Kinnaur) (Fig. 2). The Principal co-ordinate analysis (PCoA) analysis also classified the populations into four major groups, which were in line with the clustering analysis (Fig. 1b). Further, the Bayesian based STRU CTU RE prediction revealed that maximum likelihood clustering was retrieved when individuals were clustered into two (K = 2) and four groups (K = 4). The K = 2 groups revealed that, D. hatagirea populations from Kinnaur (Sangla and Chittkul), Kullu (Tosh), and Lahaul valley (Shakoli) were grouped together in one group and populations from Spiti valley (Mane, Sichling, Shego, Giyu and Rangrik) were grouped together. Further, the maximum likelihood clustering with K = 4 revealed that Bathad population from Kullu region and Rangrik population from Spiti valley belongs to different groups, however, one individual from Mane and Giyu revealed intermixing (Fig. 3).

Discussion
The drastic changes in the environmental conditions due to climate change and overexploitation by human's intervention, the risks of extinction of various endangered plants of Himalayan region has increased many times 19 . Orchids are the largest plant families which are commercially traded for various purposes in local, national and international markets. D. hatagirea is among the widely utilised and traded Orchid around the Table 3. Genetic differentiation coefficient Fst (below diagonal) and gene flow Nm (above diagonal) between 10 populations of Dactylorhiza hatagirea.  Unfortunately, each kilos of the D. hatagirea comprised of around 100 individuals and much of the procurement, harvesting and trade is done illegally and unsustainably, results in its critically endangered status with huge conservation concerns 7 . Understanding the genetic diversity of endangered plant species is of utmost importance, because high genetic diversity is crucial for maintaining the fitness and survival of natural populations, which also ensure the presence of pre-adapted genotypes 21 . Moreover, information related to genetic diversity present within and among populations of endangered plants like D. hatagirea is important to establish the conservation strategies 22 . www.nature.com/scientificreports/ In current study, genetic diversity of D. hatagirea populations from three cold dessert regions of Himachal Pradesh was investigated, because the plant is considered as model to study orchids and their immense pharmacological and therapeutic properties like anticancer, antioxidants and activation of immune system to fight infectious diseases 23 . Despite the medicinal and economical importance of D. hatagirea along with its critically endangered status, only few molecular studies using ISSR and RAPD markers have been carried out to understand the genetic diversity and population structure pattern of the D. hatagirea populations. However, SSR markers which are widely utilised and recommended for diversity analysis due to their co-dominant and high polymorphic nature has not been used till date to analyse the diversity pattern of D. hatagirea populations in the cold desert regions of North-Western Himalayas 24 .

Genetic diversity analysis based on SSR markers in D. hatagirea.
In current study, overall low levels of genetic diversity (H o = 0.185 and H e = 0.158) was recorded in the populations of D. hatagirea which was also reported in earlier studies where genetic diversity of D. hatagirea from Ladakh region was estimated using ISSR and RAPD markers 12,13 . The low genetic diversity in few populations of cold desert regions of Kullu, Kinnaur and Lahaul Spiti valley stipulates that the species is under serious threat of extinction in these regions 25 . Moreover, the low genetic diversity recorded might be because of the isolated populations, habitat specificity, www.nature.com/scientificreports/ restricted distribution and high anthropogenic activities 26 . Moreover, the populations of D. hatagirea in alpine meadows ecosystem, prefers humid laden undulating habitat and its population density in alpine meadows usually ambit between 0.70 and 2.43 individuals/m 28,26 . Another possible reason for the low genetic diversity is the small size of the populations and large spatial distance between the populations which hinders the gene flow between the populations, results in self-pollination and inbreeding 24 In case of D. hatagirea, selfing and crossing among related individuals has been reported 12 . All these factors might have resulted into the reduced life span, and health of D. hatagirea individuals with exhausted capacity to grow population size 27 . Current study has reported for the first time the genetic diversity pattern of the 10 populations representing forty-eight individuals from three cold desert regions of Himachal Pradesh. The study reveals the overall variation in the polymorphism level between these populations. Overall, genetic diversity (I, H o , H e , and PIC) at population level was recorded uniform for majority of populations and was moderated to high in only three populations viz., Mane, Rangrik, Sichling from Spiti valley and one Tosh from Kullu valley. The genetic diversity pattern in these populations indicates low level of human interferences/disturbance in comparison to other populations. Further, individuals from Spiti valley represents the hub for genetic diversity of D. hatagirea. However, low levels of genetic diversity in other populations might be due to overexploitation, habitat destruction, poor regeneration due to its specific requirement of growth association with mycorrhiza fungi for germination. Gene flow, and genetic differentiation. Two important parameters viz., the genetic differentiation coefficient and gene flow plays an important role to understand the genetic structure of the populations and these two parameters are negatively correlated 24,28 . In current study, gene flow between the population of D. hatagirea was recorded < 1 (mean N m = 0.363), which indicates that genetic differentiation is due to genetic drift 29 and it might be an important factor which is influencing the genetic structure of D. hatagirea populations in these regions of Himachal Pradesh. Moreover, Mantel test analysis revealed weak positive correlation between the genetic diversity and geographic distance between populations of D. hatagirea. Therefore, the differences in the genetic differentiation might be due to few geographic barriers, resulting in isolating populations in different gene pools of D. hatagirea. Further, AMOVA analysis revealed that high molecular variance was recorded among populations and low within populations, which further stipulates habitat fragmentation due to human interventions and anthropogenic activities, leading to low genetic diversity within populations 13 . Moreover, lower genetic variation within populations might be due to vegetative propagation of D. hatagirea in the local regions of Western Himalayas.
Genetic relationship and population structure. The PCoA, UPGMA and STRU CTU RE analysis was also performed in the study to predict the genetic diversity and population structure of 48 individuals representing 10 populations of D. hatagirea. The distribution pattern of majority individuals of different populations on the basis of these analysis revealed consistency with their geographic origin. The UPGMA analysis grouped the populations into 4 major clusters wherein, cluster I and IV represented population from Spiti valley; cluster II represented population from Kullu and cluster III revealed grouping of populations form Kinnaur, Kullu and Lahaul valley. Despite the long distance between the populations which were represented in cluster III, stipulates low levels of gene flow (mean N m = 0.363) among them, might be due to collection of D. hatagirea from one region and its introduction to the other region. Interestingly, two populations viz., Bathad and Rangrik might be of completely different origin as compare to other populations, also possessed large number of private alleles, and might have role in response to selection and of high evolutionary significance 18 . Structure analysis also clearly suggests that individuals of these populations belong to four major clusters with admixture of few individuals in which populations from Spiti valley seems to be different from rest of the populations. Furthermore, PCoA was also found similar to STRU CTU RE results which validate the UPGMA tree. Small admixture retrieved in the structure stipulates the low gene flow (mean N m = 0.363) among some individuals of 10 populations. Overall, above analysis revealed that it is most likely that majority of the D. hatagirea populations in these regions might have been distributed via anthropogenic method and propagated vegetatively, resulted in low population differentiation, gene flow and are genetically similar (Fig. 4).  (Table 6 and Fig. 4).

DNA Isolation and PCR amplification.
Total genomic DNA extraction was performed using CTAB protocol with little modifications 30 . The quantity and quality of the extracted DNA was estimated using NanoDrop 2000 OD 260 /OD 280 (Thermo Scientific, Lithuania) and integrity using 1% agarose gel with uncut λ as reference.   www.nature.com/scientificreports/ PAGE (Poly Acrylamide Gel Electrophoresis) gel and were visualised using silver staining method along with size estimation of amplified product against 50 bp ladder (ThermoFisher).
SSR data scoring, genetic diversity and population structure prediction. Stable and highly polymorphic fifteen SSR markers (KSSR) were utilised for genotyping of 48 individuals from 10 wild natural populations belonging to three cold desert regions (Kullu, Kinnaur and Lahaul & Spiti) of Himachal Pradesh (Table 7). Presence (1) and absence (0) of bands of a particular size (molecular weight) in the genotyping profiles of all the individuals were scored manually in the binary format. To establish the genetic relationship among collected individuals, each band on the genotyping profile was considered as loci or marker. Scored data for all individuals was utilised for Nei's gene diversity (h), shannon's information index (I), effective number of alleles (N e ), percent polymorphic marker loci (%PPL), analysis of molecular variance (AMOVA), principal coordinate analysis (PCoA), and mantel test of geographic and genetic distance estimation using GenAlEx version 6.5 31  . Cluster analysis was performed using Nei's genetic similarity matrix generated in GenAlex from 15 SSR markers using UPGMA (Unweighted Pair Group Method with Arithmetic mean). The data was further utilised to develop a dendrogram using unpaired group mean arithmetic and Jaccard's similarity coefficients in DARWIN 6.0 software 35 . STRU CTU RE analysis based on Bayesian clustering was performed using STRU CTU RE version 2.3.4 36 . The population genetic structure (k) was classified by six iterations and value of k was set between "1 to 5" with a burn-in phase 100,000, also Markov Chain Monte Carlo (MCMC) repeats was set 100,000 after burn in and final structure determination was performed using Structure Harvester version 0.6.94 37 .
Conservation implications. D. hatagirea is one of the most important species of Orchidaceae family which is exploited for its bioactive compounds viz., Dactylorhins A-E, and Dactyloses A-B. The species hold immense significance in the clinical and drug discovery research to identify potential bioactive compounds for treatment of diseases associated with respiratory, circulatory and reproductive systems. Unfortunately, due to www.nature.com/scientificreports/ over-exploitation and unsustainable management, the populations of D. hatagirea are dwindling in their natural habitat and loss of genetic diversity within/among these populations further increases the risk of extinction due to factors like inbreeding depression. Moreover, small size of the population also accumulates deleterious genetic variations/mutations due to genetic drift and results in highly vulnerable populations 38,39 . Therefore, increasing the population size is one of the important parameter in conservation of endangered plant species, as effective population size (N e ) have huge evolutionary potential and essential to prevent inbreeding depression 40 . The population structure and genetic diversity analysis of D. hatagirea revealed low gene flow. Interestingly, populations from two regions viz., Bathad from Kullu and Rangrik from Spiti valley revealed unique and different genetic structure, can be considered as potential populations for conservation strategy planning via introducing them to different regions to induce artificial gene flow 41 . Therefore, conservation of diverse individuals of each population is important to preserve the genetic diversity of D. hatagirea in the cold desert regions of Himachal Pradesh as diverse gene pools improves the chances of adaptation for future climatic conditions 42 .